################################################################################
# Libraries and Import
################################################################################
rm(list=ls())

# Libraries
library(sandwich)
library(tidyverse)
library(stargazer)

# Import
load(file="data/dat_analysis_JP.RData")
jp <- vignette
load(file="data/dat_analysis_CN.RData")
cn <- dat_vignette

# Set seed
set.seed(30)

################################################################################
# Clean covariates
################################################################################
# China
cn$age <- as.numeric(as.character(cn$age_group))
cn$international_sys <- as.numeric(as.character(cn$international_sys))
cn$perception_fp <- as.numeric(as.character(cn$perception_fp))
cn$critique <- as.numeric(as.character(cn$critique))
cn$China_citizen <- as.numeric(as.character(cn$China_citizen))
cn$territory <- as.numeric(as.character(cn$territory))

cn$nationalism <- 0
cn$nationalism[cn$vignette_group == 1| cn$vignette_group == 2] <- 1
cn$repression <- 0
cn$repression[cn$vignette_group == 2| cn$vignette_group == 4] <- 1

# Japan
jp$age <- as.numeric(as.character(jp$age_group))
jp$international_sys <- as.numeric(as.character(jp$international_sys))
jp$perception_fp <- as.numeric(as.character(jp$perception_fp))
jp$critique <- as.numeric(as.character(jp$critique))
jp$Japan_citizen <- as.numeric(as.character(jp$Japan_citizen))
jp$territory <- as.numeric(as.character(jp$territory))

jp$nationalism <- 0
jp$nationalism[jp$vignette_group == 1| jp$vignette_group == 2] <- 1
jp$repression <- 0
jp$repression[jp$vignette_group == 2| jp$vignette_group == 4] <- 1

################################################################################
# Split each country sample into two datasets by treatment group
################################################################################
# Split dataset by treatment group - Japan
nat_jp = jp %>% filter(!is.na(nat)) %>% rename(z = nat)
econ_jp = jp %>% filter(!is.na(econ)) %>% rename(z = econ)

# Convert China response variable from factor to numeric
cn$D1 = as.numeric(as.character(cn$D1))

# Split dataset by treatment group - China
nat_cn = cn %>% filter(vignette_group == 1 | vignette_group == 2)
econ_cn = cn %>% filter(vignette_group == 3 | vignette_group == 4)

# Define treatment variable as "z"
nat_cn$z = with(nat_cn, ifelse(vignette_group == 1, 0, 1))
econ_cn$z = with(econ_cn, ifelse(vignette_group == 3, 0, 1))

###############################################################################
# Estimation of heterogenous effects: China/Japan by Party
###############################################################################
# Define ruling party in Japan as LDP/Komeito coalition
nat_jp$party <- with(nat_jp, ifelse((party == 5 | party == 9), 1, 0))
econ_jp$party <- with(econ_jp, ifelse((party == 5 | party == 9), 1, 0))

# Create baseline tables: estimate models with no covariates
nat_cn_party = lm(D1 ~ z + party + z:party, data = nat_cn)
nat_jp_party = lm(D1 ~ z + party + z:party, data = nat_jp)
econ_cn_party = lm(D1 ~ z + party + z:party, data = econ_cn)
econ_jp_party = lm(D1 ~ z + party + z:party, data = econ_jp)

# Create baseline tables: estimate models with covariates
nat_cn_party_cov = lm(D1 ~ z + party + z:party + 
                           age + past_resi + current_resi, 
                  data = nat_cn)

nat_jp_party_cov = lm(D1 ~ z + party + z:party 
                             + male + territory, 
                      data = nat_jp)

summary(nat_jp_party_cov)

econ_cn_party_cov = lm(D1 ~ z + party + z:party + 
                         age_group + regime_insider + past_resi + current_resi, 
                  data = econ_cn)

econ_jp_party_cov = lm(D1 ~ z + party + z:party + age_group + male + income, 
                       data = econ_jp)

###############################################################################
# Tables of heterogenous effects: China/Japan by Party
###############################################################################
# Calculate robust standard errors for use in table output
nat_cn_party_se <- sqrt(diag(vcovHC(nat_cn_party, type = "HC")))
econ_cn_party_se <- sqrt(diag(vcovHC(econ_cn_party, type = "HC")))
nat_jp_party_se <- sqrt(diag(vcovHC(nat_jp_party, type = "HC")))
econ_jp_party_se <- sqrt(diag(vcovHC(econ_jp_party, type = "HC")))

# Create stargazer output
stargazer(nat_cn_party, econ_cn_party, nat_jp_party, econ_jp_party,
          se = list(nat_cn_party_se, econ_cn_party_se, nat_jp_party_se, econ_jp_party_se),
          #p = makerobustpslist(nat_cn_party, econ_cn_party, nat_jp_party, econ_jp_party), 
          dep.var.labels = "Government support",
          covariate.labels = c("Government interence",
                               "Ruling party affiliate",
                               "Interference x ruling party affiliate",
                               "Constant"),
          column.labels=c("China: Nationalist", "China: Economic", 
                          "Japan: Nationalist", "Japan: Economic"),
          title = "The Baseline Result",
          label = "tab:nat_econ_party",
          align = TRUE,
          header = TRUE,
          omit.stat = c("f", "rsq", "ser"), 
          out = "tables/T3. vignette_party.tex")

###############################################################################
# Charts of heterogenous effects: China/Japan by Party
###############################################################################
# Extract ATE and SEs from regression models
nat_jp_ldp = lm(D1 ~ z, data = nat_jp, subset = party == 1)
nat_jp_non_ldp = lm(D1 ~ z, data = nat_jp, subset = party == 0)
nat_cn_ccp = lm(D1 ~ z, data = nat_cn, subset = party == 1)
nat_cn_non_ccp = lm(D1 ~ z, data = nat_cn, subset = party == 0)
econ_jp_ldp = lm(D1 ~ z, data = econ_jp, subset = party == 1)
econ_jp_non_ldp = lm(D1 ~ z, data = econ_jp, subset = party == 0)
econ_cn_ccp = lm(D1 ~ z, data = econ_cn, subset = party == 1)
econ_cn_non_ccp = lm(D1 ~ z, data = econ_cn, subset = party == 0)

ate <- vector(mode = "numeric", length = 8)
se <- vector(mode = "numeric", length = 8)
party <- vector(mode = "character", length = 8)
country <- vector(mode = "character", length = 8)

ate[1] = summary(nat_jp_ldp)$coefficients[2]
ate[2] = summary(nat_jp_non_ldp)$coefficients[2]
ate[3] = summary(nat_cn_ccp)$coefficients[2]
ate[4] = summary(nat_cn_non_ccp)$coefficients[2]
ate[5] = summary(econ_jp_ldp)$coefficients[2]
ate[6] = summary(econ_jp_non_ldp)$coefficients[2]
ate[7] = summary(econ_cn_ccp)$coefficients[2]
ate[8] = summary(econ_cn_non_ccp)$coefficients[2]

se[1] <- sqrt(diag(vcovHC(nat_jp_ldp, type = "HC")))[2]
se[2] <- sqrt(diag(vcovHC(nat_jp_non_ldp, type = "HC")))[2]
se[3] <- sqrt(diag(vcovHC(nat_cn_ccp, type = "HC")))[2]
se[4] <- sqrt(diag(vcovHC(nat_cn_non_ccp, type = "HC")))[2]
se[5] <- sqrt(diag(vcovHC(econ_jp_ldp, type = "HC")))[2]
se[6] <- sqrt(diag(vcovHC(econ_jp_non_ldp, type = "HC")))[2]
se[7] <- sqrt(diag(vcovHC(econ_cn_ccp, type = "HC")))[2]
se[8] <- sqrt(diag(vcovHC(econ_cn_non_ccp, type = "HC")))[2]

party[1] = "Ruling Party"
party[2] = "Non Ruling Party"
party[3] = "Ruling Party"
party[4] = "Non Ruling Party"
party[5] = "Ruling Party"
party[6] = "Non Ruling Party"
party[7] = "Ruling Party"
party[8] = "Non Ruling Party"

country[1] = "Japan: Nationalist"
country[2] = "Japan: Nationalist"
country[3] = "China: Nationalist"
country[4] = "China: Nationalist"
country[5] = "Japan: Economic"
country[6] = "Japan: Economic"
country[7] = "China: Economic"
country[8] = "China: Economic"

party_ate = data.frame(ate, se, party, country)

# Create chart
ate <- ggplot(subset(party_ate, 
              country %in% c("Japan: Nationalist" , "China: Nationalist")), 
  aes(party, ate, color = party)) +
  geom_point(size = 3, alpha = 0.9) + 
  scale_color_manual(values = c("steelblue2", "seagreen3"), guide = FALSE) +
  geom_errorbar(aes(ymin = ate - 1.96*se, ymax = ate + 1.96*se),
                position = position_dodge(),
                color="black", alpha = 0.3,
                size=0.5,
                width = 0.3) +
  geom_text(aes(label = round(ate, 2)), nudge_x = 0.25, 
            size = 3, color = "black") +
  facet_wrap( ~ country, ncol = 2) + 
  geom_hline(yintercept = 0, linetype = "dashed") +
  xlab("Party Affiliation") +
  ylab("Mean change in government approval") + 
  theme_classic() +
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(axis.text=element_text(size = 10)) +
  theme(legend.position = "bottom")

ggsave("figs/7. ate_china_japan_party_nationalism.pdf", ate, height = 5, width = 6)

ate <- ggplot(subset(party_ate, country %in% c("Japan: Economic" , "China: Economic")), 
              aes(party, ate, color = party)) +
  geom_point(size = 3, alpha = 0.9) + 
  scale_color_manual(values = c("steelblue2", "seagreen3"), guide = FALSE) +
  geom_errorbar(aes(ymin = ate - 1.96*se, ymax = ate + 1.96*se),
                position = position_dodge(),
                color="black", alpha = 0.3,
                size=0.5,
                width = 0.3) +
  geom_text(aes(label = round(ate, 2)), nudge_x = 0.25, 
            size = 3, color = "black") +
  facet_wrap( ~ country, ncol = 2) + 
  geom_hline(yintercept = 0, linetype = "dashed") +
  scale_y_continuous(limits = c(-1, 0.05)) +
  xlab("Party Affiliation") +
  ylab("Mean change in government approval") + 
  theme_classic() +
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(axis.text=element_text(size = 10)) +
  theme(legend.position = "bottom")

ggsave("figs/8. ate_china_japan_party_economic.pdf", 
       ate, height = 5, width = 6)